/*
 * MatAlgebra.cpp
 *
 *  Created on: May 27, 2011
 *      Author: nzhao
 */

#include "share/Usrdef.h"
#include "share/MatAlgebra.h"


void realsym_mat_diag(int dim, double * mat, gsl_vector * eval, gsl_matrix * evec)
{
	gsl_matrix_view m
		= gsl_matrix_view_array (mat, dim, dim);

	gsl_eigen_symmv_workspace * w =	gsl_eigen_symmv_alloc (dim);

	gsl_eigen_symmv (&m.matrix, eval, evec, w);
	gsl_eigen_symmv_free (w);
	gsl_eigen_symmv_sort (eval, evec,GSL_EIGEN_SORT_VAL_ASC);

}

void hermite_mat_diag(int dim, double * mat, gsl_vector * eval, gsl_matrix_complex * evec)
{
	gsl_matrix_complex_view m
		= gsl_matrix_complex_view_array (mat, dim, dim);

	gsl_eigen_hermv_workspace * w =	gsl_eigen_hermv_alloc (dim);

	gsl_eigen_hermv (&m.matrix, eval, evec, w);
	gsl_eigen_hermv_free (w);
	gsl_eigen_hermv_sort (eval, evec,GSL_EIGEN_SORT_VAL_ASC);

}
void MKL_mat_mult(int dim, char * transa, char * transb, CLD alpha, CLD * mat1, CLD * mat2, CLD beta, CLD * mat_res)
{
    int  m, n, k, lda, ldb, ldc, i, j;


    m = dim;
    n = dim;
    k = dim;
    lda = dim;
    ldb = dim;
    ldc = dim;

    zgemm(transa, transb, &m, &n, &k, &alpha, mat1, &lda, mat2, &ldb, &beta, mat_res, &ldc);
//    for(i = 0; i < dim; i ++)
//    {
//    	for(j = 0; j < dim; j ++)
//    	{
//    		cout << mat_res[i*dim + j] << "\t";
//    	}
//    	cout << endl;
//    }

}

void MKL_real_mat_mult(int dim, char * transa, char * transb, double alpha, double * mat1, double * mat2, double beta, double * mat_res)
{
    int  m, n, k, lda, ldb, ldc, i, j;


    m = dim;
    n = dim;
    k = dim;
    lda = dim;
    ldb = dim;
    ldc = dim;

    dgemm(transa, transb, &m, &n, &k, &alpha, mat1, &lda, mat2, &ldb, &beta, mat_res, &ldc);
}

void MKL_mat_vec_mult(int dim, char * trans,CLD alpha, CLD * mat, CLD * vect, CLD beta, CLD * vect_res)
{
    int  m, n, lda, incx, incy, i, j;


    m = dim;
    n = dim;
    lda = dim;
    incx = 1;
    incy = 1;

    zgemv(trans, &m, &n, &alpha, mat, &lda, vect, &incx, &beta, vect_res, &incy);
}

void MKL_matrix_copy(char trans, int dim, CLD alpha, CLD * src, CLD * dst)
{
	int i, j;
    char ordering= 'C';	//Column major
    int src_stride = dim;
    int dst_stride = dim;

    mkl_zomatcopy(ordering, trans, dim, dim, alpha, src, src_stride, dst, dst_stride);

//    //src
//    cout << "src mat:" << endl;
//    for(i = 0; i < dim; i ++)
//    {
//    	for(j = 0; j < dim; j ++)
//    	{
//    		cout << src[i*dim + j] << "\t";
//    	}
//    	cout << endl;
//    }
//    cout << endl;
//
//    //dst
//    cout << "dst mat:" << endl;
//    for(i = 0; i < dim; i ++)
//    {
//    	for(j = 0; j < dim; j ++)
//    	{
//    		cout << dst[i*dim + j] << "\t";
//    	}
//    	cout << endl;
//    }
//    cout << endl;
}


void MKL_hermite_mat_diag(int dim, CLD * mat, double * eval, CLD * evec)
{
   /* Locals */
    int n = dim, lda = dim, ldz = dim, il, iu, m, info, lwork, lrwork, liwork;
    double abstol, vl, vu;
    int iwkopt;
    int* iwork;
    double rwkopt;
    double* rwork;
    CLD wkopt;
    CLD * work;

    char jobz[]  = "Vectors";
    char range[] = "All";
    char uplo[]  = "Lower";

    /* Local arrays */
    int * isuppz = new int [2*dim];
    double * w = eval;
    CLD * z = evec;
    CLD * a = mat;

    /* Negative abstol means using the default value */
    abstol = -1.0;
    /* Query and allocate the optimal workspace */
    lwork = -1;
    lrwork = -1;
    liwork = -1;

//    cout << "querry..." << endl;
    zheevr( jobz, range, uplo, &n, a, &lda, &vl, &vu, &il, &iu,
                        &abstol, &m, w, z, &ldz, isuppz, &wkopt, &lwork, &rwkopt, &lrwork,
                        &iwkopt, &liwork, &info );
    lwork = (int) wkopt.real();
    work = new CLD [lwork];
    lrwork = (int)rwkopt;
    rwork = new double [lrwork];
    liwork = iwkopt;
    iwork = new int [liwork];

    /* Solve eigenproblem */
    zheevr( jobz, range, uplo, &n, a, &lda, &vl, &vu, &il, &iu,
                    &abstol, &m, w, z, &ldz, isuppz, work, &lwork, rwork, &lrwork,
                    iwork, &liwork, &info );

//    int i, j;
//    for(i = 0; i < n; i++)
//    {
//        cout << w[i] << "\t {";
//        for(j = 0; j < n; j++)
//            cout << z[ i * n + j ] << ",";
//        cout << "}" << endl;
//    }
//    cout << endl;

    delete [] work;
    delete [] rwork;
    delete [] iwork;
}

void MKL_matrix_complex_row_transform(int dim, CLD * exp_vals, CLD * right_mat)
{
	int i, j;
	for(i = 0; i < dim; i ++)
	{
		for(j = 0; j < dim; j ++)
		{
			right_mat[j*dim + i] *= exp_vals[i];
		}
	}

}

void MKL_i_t_Hermite_Matrix_Exp(int dim, double * eigvals, CLD * eigvec, double t, CLD * res_mat)
{
	int i, j;
	CLD * i_eigen_val 		= new CLD [dim];
	CLD * exp_i_eigen_val 	= new CLD [dim];
	CLD * res_mat1 			= new CLD [dim * dim];

	for(i = 0; i < dim; i++)
		i_eigen_val[i] = CLD (0.0, t * eigvals[i]);
	vzExp(dim,i_eigen_val,exp_i_eigen_val);

	MKL_matrix_copy('C', dim, COMPLEX_ONE, eigvec, res_mat1);
	MKL_matrix_complex_row_transform(dim, exp_i_eigen_val, res_mat1);
	MKL_mat_mult(dim, "NoTranspose", "NoTranspose", COMPLEX_ONE, eigvec, res_mat1, COMPLEX_ZERO, res_mat);

	delete [] i_eigen_val;
	delete [] exp_i_eigen_val;
	delete [] res_mat1;
}
